Bayesian Analysis of QENS data: Prom parameter determination to model selection 
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The extraction of any physical information from quasielastic neutron scattering spectra is generally 
done by fitting a model to the data by means of x 2 minimization procedure. However, as pointed 
out by the pioneering work of D.S. Sivia et al. [J|, also another probabilistic approach based on 
Bayes theorem |2| can be employed. In a nutshell, the main difference between the classical \ 2 
minimization and the Bayesian approach is the way of expressing the final results: In the first case, 
the result is a set of values of parameters with a symmetric error [Pi ± Si) and a figure of merit 
such as x 2 1 whereas in the second case the results are presented as probability distribution functions 
(PDF) of both, parameters and merit figure. In this contribution, we demonstrate how final PDFs 
are obtained by exploring all possible combinations of parameters that are compatible with the 
experimental error. This is achieved by allowing the fitting procedure to wander in the parameter 
space with a probability of visiting a certain point P = exp(— x 2 the so called Gibbs sampling. 
Three advantages of this method will be emphasized: First, correlations between parameters are 
automatically taken into account, which implies, for example, that parameter errors are correctly 
calculated, correlations show up in a natural way and ill defined parameters (i.e. parameters for 
which data only support the calculation of lower or upper bounds) are immediately recognized from 
their PDF. Second, it is possible to calculate the likelihood of a determined physical model, and 
therefore to select the one among many that fits the data best with a minimal number of parameters, 
in a correctly defined probabilistic way. Finally, in the case of a low count rate where the Gaussian 
approximation to the Poisson statistics fails, this method can also be used by simply redefining x 2 ■ 



INTRODUCTION 

Science is based on the success of an hypothesis to de- 
scribe experimental results, i.e. is based on the amount 
of " truth" and " falsity" of an hypothesis when contrasted 
with experimental results [3[. In order to find a quanti- 
tative method to determine this " amount of truth" , hy- 
potheses in science should at the end be reduced to a 
mathematical expression depending on a set of parame- 
ters with some physical meaning. The "amount of truth" 
is then determined by fitting the mathematical model to 
some experimental data. The general method to do so 
is to minimize the squared distance between experimen- 
tal data and the points generated by the mathematical 
model. Furthermore, taking also into account the error 
associated with experimental data, a figure of merit x 2 
can be defined 
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where n is the number of experimental points and m 
is the number of parameters, D k (fc=l,... n) are the 
experimental data, H^{Pi} (k=l,... n) are the values 
obtained from our hypothesis (the mathematical model) 
using the {Pi} (i=l,... m) set of parameters contained 
in the model, and <7fc (fc=l,... n) are the experimental 
errors associated with the respective Dk- Therefore, the 
fitting procedure has a twofold goal: first, to find the 
set of parameters {Pi} which describes the experimen- 
tal data best, and second, using this set of parameters, 



to define a figure of merit which quantifies the "amount 
of truth" of the proposed hypothesis. In order to be 
able to compare different hypotheses with different num- 
bers of parameters it is reasonable to define a figure of 
merit which penalizes additional parameters such as the 
reduced \ 2 defined as: 
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where n is the number of experimental points and m is 
the number of parameters, so n — m is the number of 
degrees of freedom. The aforementioned way to quantify 
how good experimental data are described by a hypoth- 
esis is based on what is called a " frequentist" approxi- 
mation to the problem and has many drawbacks as- 
sociated with both the fitting procedure and the way to 
quantify the correctness of the hypothesis describing ex- 
perimental data. 

Data fitting is usually done by minimizing the 
aforementioned % 2 (equation [T]) using the Levenberg- 
Marquardt algorithm, which aims to find the minimum 
of the x 2 {Pi} hypersurface. Unfortunately, usually local 
minima make the algorithm unable to find the absolute 
minimum. For this reason, this method can find a final 
solution only when the algorithm is initialized with pa- 
rameters near the global minimum. The final solution is 
then characterized by a set of parameters with an associ- 
ated error (Pk±£k) and the figure of merit xt- This way 
of quantifying the best fit to the data is based on the sup- 
position that there is only one minimum in the x 2 (Pk) 



2 



hypersurface compatilbe with data error, and that the 
functional dependence of x 2 {Pk) is quadratic on each pa- 
rameter (i. e. one can stop at the second term of a Taylor 
expansion of the obtained minimum), and thus allowing 
only symmetric errors. Moreover, errors are usually cal- 
culated disregarding possible correlations between them 
Q and are thus generally underestimated. 

We present in this work a method both to perform 
fittings and to analyze results based exclusively on prob- 
ability by using what is called Baycsian inference. The 
main difference with the previously exposed frequentist 
method is the absence of any supposition on the x 2 {Pi} 
landscape which will rather be explored using the proba- 
bilities determined from experimental data. The method 
results in a different way to express fitted parameters 
and the figure of merit showing all the complexity of 
the final solution: they become Probability Distribution 
Functions (PDFs) obtained directly from exploring the 
X 2 {Pk} hypersurface. 

The paper will be organized as follows: first the ubiq- 
uitous x 2 w iU be defined using exclusively probability 
theory, and on this basis a method to sample the X 2 {Pk] 
hypersurface will be presented: the Gibbs sampling. We 
will then refer on how both the frequentist and Bayesian 
methods select an hypotheses among others, stressing the 
advantages of using the second approach. Finally the 
presented method implemented in the FABADA package 
[|| will be applied to three real cases related to neutron 
scattering each stressing different aspects of the proposed 
method. In the first example, the importance of let- 
ting parameters free or fixed in the fitting process will 
be stressed. The second example will focus on the PDF 
obtained from a set of data fitted simultaneously, and 
model selection will be addressed in the third example. 



WHAT IS BEHIND THE UBIQUITOUS x 2 ? 

The objective of the so called Bayesian methods [l], Q 
is to find the probability that a hypothesis is true given 
some experimental evidence. This is done taking into ac- 
count both our prior state of knowledge concerning the 
hypothesis, and the likelihood that the data is described 
by the proposed hypothesis. Using probability notation, 
and only considering the case that the experiment con- 
sists of a series of data D k and that the hypothesis is rep- 
resented by H k , we can relate the aforementioned prob- 
abilities using the Bayes theorem Q: 
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where P(H k \ D k ) is called the posterior, the proba- 
bility that the hypothesis is in fact describing the data. 
P{D k | H k ) is named the likelihood, the probability that 
our data is well described by our hypothesis. P(H k ) is 
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FIG. 1: Poisson statistics followed by a counting experiment 
such as a neutron scattering experiment (lines). For an in- 
creasing number of counts, the Poisson distribution can be 
approximated by a Gaussian function (points) with rr = y/n, 
being n the number of counts. 



called the prior, the knowledge we have beforehand about 
the hypothesis, and P(D k ) is a normalization factor to 
assure that the integrated posterior probability is unity. 
In the method here presented we will assume no prior 
knowledge (maximum ignorance prior [2j), and in this 
special case Bayes theorem takes the simple form: 
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where L is a short notation for likelihood. 

We need first to find the likelihood that one data point 
D k is described by the mathematically modeled hypoth- 
esis H k . In a counting experiment such as those related 
to neutron scattering this probability follows a Poisson 
distribution 
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Nevertheless, for a high enough number of counts, the 
Poisson PDF can be well approximated by a Gaussian 
one with a = V ' D k as it is shown in figure [T] and hence 
the likelihood that the set of data points D k is correctly 
described by the hypothesis H k can be written as 
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Therefore, we have found the meaning of the ubiqui- 
tous x 2 based only on probabilistic grounds: it is related 
to the probability that a certain set of data is well de- 
scribed by an hypothesis, and hence the goal of mini- 
mizing x 2 is finding a set of parameters that maximizes 
the likelihood associated with the proposed mathematical 
model. The probability theory behind x 2 allows there- 
fore also to deal with the case of experiments with only 
few counts where the Gaussian approximation is not valid 
anymore and the Poisson distribution must be employed 
simply by redefining x 2 as 
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THE BAYESIAN METHOD: GIBBS SAMPLING 
OF PARAMETER SPACE 

The probabilistic understanding of % 2 makes it possible 
to define a unique method, first to fit the experimental 
data, and then to analyze the obtained results, using a 
Markov Chain Monte Carlo (MCMC) technique. A set 
of parameters P™ ew is generated from an old set P° ld 
by randomly changing one of the parameters The 
probability to accept the new set of parameters is given 
by 
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where Xnew an( l Xoid correspond to the % 2 (as defined 
in equation [T]) for the new and old set of parameters. 
This way of exploring the parameter space (called Gibbs 
sampling) is similar to the way used to find the pos- 
sible molecular configurations of a determined system 
at a given temperature using the classical Monte Carlo 
method: the values of physical constants such as the po- 
tential energy will in fact be a PDF related to all the 
configurations explored by the Monte Carlo method. It 
is therefore possible to relate energy to J2(^k~Dk) 2 , the 
magnitude giving information about the fit quality of the 
hypothesis with respect to the data, and temperature to 
the error associated with the data (T ~ a 2 ). This wat of 
exploring the parameter space has two main advantages: 

• In the fitting process, the Bayesian method is able 
to accept a new set of parameters that do not de- 
crease x 2 j if this change is compatible with the ex- 
perimental error and therefore does not get stuck 
in local minima as the Levenberg-Marquardt algo- 
rithm. In other words, the presented method is 
able to go "uphill" in the X 2 {Pi} hypersurface if 
the barrier is compatible with the error. Neverthe- 
less, in order to avoid the presented algorithm to 
get stuck even in the case when barriers are greater 



than those associated to experimental error a sim- 
ulated annealing can be used. This algorithm cal- 
culates a fictitious x 2 = 2/s=i wriere T 
is a constant defined to artificially increase the ex- 
perimental error, and by similitude with classical 
Montecarlo simulation is named as " temperature" . 
Fittings are then started at high temperature, and 
the system is relaxed by lowering the temperature 
up to T=l. 

• Concerning the analysis of the results obtained by 
the fitting, the exploration of the whole parame- 
ter space compatible with data using the MCMC 
method allows both to find the PDF associated 
with the likelihood directly related the figure of 
merit x 2 ( see equation |4|), and the parameters, 
taking into account possible correlations between 
them, or minima not describable by a quadratic 
approximation. 



MODEL SELECTION 

Data can usually be described by more than one hy- 
pothesis, each implying a different physical mechanism 
to explain experimental results. Albeit the importance 
to perform model selection accurately, vague arguments 
are usually given to prefer a model among others and 
usually no quantitative arguments are given to justify 
why an hypothesis is preferred, although it is possible to 
do so using both the frequentist and Bayesian methods. 
Model selection can be performed using the frequentist 
approach by using the xl figure of merit (see equation HJ 
which takes into account the addition of parameters to a 
model by dividing x 2 by the degrees of freedom. There- 
fore, if two models fit the data with equal success, i.e. 
with the same % 2 , the model with less parameters (with 
the smallest xV) w iH be favored. In some sense this is 
nothing but quantifying the Ockham's razor principle: 
it is necessary to shave away unnecessary assumptions 
(parameters). Model selection performed by using xl 
has the same drawbacks as the determination of param- 
eter errors: we suppose that there is a single minimum 
in X 2 {Pi}i that this minimum parabolic depends on all 
parameters and that there are no correlations between 
parameters. In fact, if these three suppositions are ac- 
complished, then the PDF of the x 2 reads Q 



P(x 2 )oc( X 2 ) Ar/2 - 1 exp(- X 2 /2) 
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N is, in this simple case, the number of parameters. In 
figure [2] the chi-square distribution for increasing degrees 
of freedom (number of parameters) is shown. As can be 
seen in the inset from figure [2) this distribution has a 
term which is independent from the number of parame- 
ters, exp(— x 2 /2), and that decreases together with the 
quality of the fitting, or when the error associated with 
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FIG. 2: x 2 distribution ((x 2 )^ 2 " 1 exp(- X 7 2 )) for an in- 
creasing number of fitting parameters N. The inset shows the 
terms associated to the quality of the fit exp(— x 2 /2) together 
with the one depending on the number of parameters of the 
model (x 2 )^ 2 " 1 . 



the experimental data <7fc increases. The term {x 2 ) N/2 ^, 
depending on the number of degrees of freedom, increases 
exponentially with the number of parameters, displac- 
ing the maximum of the x 2 distribution to higher values. 
Therefore, even using the frequentist approach, the afore- 
mentioned preference for models that fit equally well the 
data with the minimum number of parameters is based 
on probability theory: those models with the maximum 
in the x 2 distribution placed at lower values will be pre- 
ferred. The Bayesian method finds in a natural way the 
PDF of x 2 by exploring the parameter space without 
the suppositions made in the frequentist approximation, 
hence the obtained PDF will in general not follow the x 2 
distribution described by equation [5J 



of the first molecular liquids studied by diffraction meth- 
ods). The determination of its molecular structure im- 
plies to obtain the distance between carbon and chlorine 
atoms, the Cl-Cl distance is fixed by the tetrahedral sym- 
metry, and the MSD between chlorine atoms and carbon 
and chlorine atoms. 

Experiments were performed at the diffractometer Dlb 
in the Institute Laue Langevin (Grenoble, France) using 
a wavelength of A = 2.52A (see 0). Figure [3^ shows 
a good agreement between experimental data and two 
fittings of equation [9j one with a fixed scale factor h, 
and the other with h as a free parameter. Figures [3)3, c 
show the PDF from parameters rcci and ^cci obtained 
through the two aforementioned fittings. Concerning the 
rcci PDF, we can immediately see that its determina- 
tion is robust since both fixed and free scale factor h 
lead to the same PDF. On the contrary, the PDF associ- 
ated with Icci is sensible to the way we have performed 
the fitting: if the scale factor is fixed we obtain a most 
probable value for this parameter (Icci = 0.066A), but 
for a free scale factor h only a maximum value for Zcci 
can be obtained due to the correlation between both pa- 
rameters (h and Icci)- Defining the upper limit as that 
for which the integrated probability is 0.682 (as errors 
are usually defined in the frequentist approach [ill ]) the 
upper limit Iqci = 0.02A can be determined from the 
cumulative distribution function (see fig. 0). 

This example shows the main difference compared to 
the frequentist approximation: the results are presented 
as PDF. This has the advantage that, as it happens with 
the determination of Icci leaving h free, the result to our 
parameter determination can be expressed as a limit for 
the parameter, which is impossible with the frequentist 
approximation. 



Parameter estimation: isotropic rotation 



EXAMPLES 

Determining the intramolecular structure of CCI4 

Molecular structure can be calculated from diffraction 
experiments by fitting the high q-range of the scattering 
function S(q) to the following equation (see Q): 



Quasi Elastic Neutron Scattering (QENS) is perfectly 
suited to determine the molecular dynamics in the liquid 
phase. Usually this dynamics is studied by splitting the 
spectra into diffusion and rotation contributions 



S(q,u>) = S(q, w) tra n S ® S(q, U>)r 



(10) 
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where bf h are the coherent cross sections for each ele- 
ment, Tij are the intramolecular distances and I 2 = (itf^) 
are the vibrational Mean Square Displacements (MSD) 
between elements i and j, and h is a scaling factor. 

In the proposed example, our objective is to calculate 
the intramolecular structure of carbon tetrachloride (one 



where S(q, w)trans is the translational contribution and 
S(q, u>)rot is associated with the rotation of the molecule, 
assuming that both movements are independent from 
each other. If we assume that the translation is described 
by a diffusion mechanism and, therefore, described by the 
Fick equation, and that rotation is isotropic [lij ]: 
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FIG. 3: (a) High momentum transfer scattering function for 
CCU, where the main contribution to S(q) is that related to 
the molecular structure. Lines are the best fits to S(q) set- 
ting the scale factor fixed and free in the fitting procedure, (b) 
Probability distribution functions obtained for the distances 
between carbon and chlorine atoms rcci, an d Mean Square 
Displacement between carbon and chlorine atom (icci) for 
both cases (fixed h, points, and free h, line). PDF(lcci) has 
been represented in logarithmic scale to show the different 
length scales explored by the Bayesian method. In addition, 
for Zcci using a h free fitting the integrated probability is 
shown, arrow points an integrated probability of 0.68 follow- 
ing the standard definition of errors in the frequentist approx- 
imation frill. 



S(q,u) Iot = A Q (q-R)S(uj) 
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where D and D r are the translational and rotational dif- 
fusion coefficients, Ai(q ■ R) are spherical bessel func- 
tions and R is the radius of rotation. We have per- 
formed QENS experiments at the TOFTOF spectrom- 
eter [l3| at the FRM II reactor (Munich) in order to de- 
termine the dynamics of 1,2-trans-dichloroethylene. The 
data were corrected for self-absorption using the FRIDA 
software 14|. A series of fittings for each temperature 
with the model described by equations [TU1 EH an d EH 
Usually, each q value is fitted separately, obtaining the 
diffusion coefficient from a second linear fit to the q 2 
dependence of the broadening of the central lorentzian 
rtrans(g 2 ) = D ■ q 2 , and the radius of rotation from 
the obtained Ao(q ■ R) = sin(qR)/qR or eventually by 
independently fitting all spectra using equation [TU] to 
each S(q = qi,uj). However, our hypothesis is described 
by the whole set of the aforementioned equations 1101 111 
and [TJl and thus errors arising from the two-step fitting 
procedure can be minimized by simply fitting the spec- 
tra S(q, oj) for all q values, i.e. fitting the complete q- 
dependent data set with only D, D r and R as physical pa- 
rameters. The results for the radius of rotation are shown 
in figure fusing the presented Bayesian method together 
with those obtained using a Levenberg-Marquardt algo- 
rithm for each q-value. First of all, because fittings were 
performed by the frequentist approximation separately 
for each q- value (see figure |4p), the radius of rotation 
has a g-dependence which is not present in the Bayesian 
fitting (see figure H^i), consequently stressing the impor- 




FIG. 4: Radius of rotation obtained (a) from the Bayesian 
fitting procedure applied to all S(q, oj) spectra and (b) from 
the frequentist fitting to each q- value individually (here only 
q — OAA and q = 0.9 A values are shown). An error bar 
obtained from the Bayesian method is plotted in (b) for S(q = 
0.4, uj) at T = WOK. In (c) the PDF for the radius of rotation 
is also shown as a function of the temperature (P(R, T)) with 
two cuts (P(R)) for T = 240K and T = 290K. 



tance of fitting the whole data set together. A fitting 
using the Bayesian algortihm has also been performed 
to a spectrum for q — OAA^ 1 and T = 300K in order 
to compare the error bars obtained by both methods. 
The Error bar using the Bayesian approach was calcu- 
lated by obtaining the PDF for the radius of rotation 
and then fitting a Gaussian function with a = e, be- 
ing e the frequentist parameter error. This error bar is 
plotted in figure |4Jd, together with that determined by 
the frequentist method. As it can be seen in the figure 
the error obtained by the presented method is much big- 
ger that that estimated by the frequentist method. The 
presented Bayesian method is therefore able to deal with 
simultaneous fitting of various curves, obtaining the PDF 
of physical parameters as a function of temperature (see 
figure Eb,d,e). 



Model selection: Diffusion in phospholipid 
membranes 



Phospholipids are the main component of cell walls 
and can also be used in technological applications as for 
example drug delivery or food industry. Their dynamics 
is studied on many time- and length-scales with different 
techniques, among them quasielastic time-of-flight neu- 
tron scattering which probes the motions that dominate 
on times of about 100 ps. 

As will be discussed in detail elsewhere [15] . the ques- 
tion arose from previous neutron scattering experiments 
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associated with x 2 is for any combination of parameters 
smaller than the one of the delta model. 

Therefore, the model comparison between the two pos- 
sibilities of broadened and non-broadened central line 
favours the model with a broadened line. 



SUMMARY 



FIG. 5: (a) Spectra of the hydrated phospholipid DMPC to- 
gether with a fit and the instrumental resolution. Internal mo- 
tions were approximated with two Lorentzians, the long-range 
motion was either assumed to be invisible (delta-shaped cen- 
tral line, not shown) or visible (broadened central line, shown 
here), (b) The \ 2 PDF associated with both the broadened 
central line and a delta function, showing that for any combi- 
nation of parameters the broadened model is preferable com- 
pared to the delta model: The data justify the assumption 
that the long-range motion is visible. 



1G, [r?j ] whether the long-range motion of phospholipids 
is visible on these times or if the motion appears rather 
localized, trapped in a cage of neighbours. This differ- 
ence can be seen in the line shape of S(q,ui): Motions 
that are localized during the observation time cause a 
central line that is not broadened beyond the resolution 
of the instrument but cause a foot in the spectrum. In 
contrast, long-range motions do broaden the central line. 

The neutron scattering experiments were performed 
with the phospholipid DMPC (1,2-Dimyristoyl-sn- 
Glycero-3-Phosphocholine) in a liquid crystal fully hy- 
drated with D2O at the neutron time-of- flight spectrom- 
eter TOFTOF at the FRM II (Munich). A typical spec- 
trum is shown in figure after standard corrections in- 
cluding self absorption and subtraction of the D2O spec- 
tra, obtained with the program FPJDA [l4[ ■ It is possi- 
ble to fit the data "satisfactorily" with both, a broadened 
and a delta-shaped central line. 

As stated before, Bayesian analysis is able to quan- 
tify how "satisfactory" the fits are, taking into account 
the whole X 2 {Pi} landscape and avoiding assumptions 
about it. In figure [SJd, the PDFs associated with x 2 f° r 
the two models are displayed. The normal Levenberg- 
Marquardt algorithm would simply return the parame- 
ters at the minimal reachable value of \ 2 together with 
this quantity. It is obvious that introducing an additional 
parameter, the nonzero width of the central line, reduces 
the x 2 ■ The question that needs to be answered is if this 
reduction is significant enough to justify the additional 
parameter. 

The \ 2 gives this answer, however relying on the as- 
sumptions discussed above. Employing the Bayesian 
Analysis, no assumptions are made as the x 2 {Pi\ land- 
scape is rendered explicitly. One can see in figure Eb that 
the model incorporating a broadened central line does 
not only yield the smaller x 2 minimum but also the PDF 



We have proposed a general Bayesian method to fit 
data, analyze results from the fit, and from these re- 
sults to perform model selection between competing hy- 
potheses. In contrast to the classical frequentist ap- 
proach, where some assumptions are done concerning the 
X 2 landscape (there is only a minimum of X 2 {-Pi} able to 
describe data within its error, this minimum has a square 
dependence on the parameters, and parameters are not 
correlated) , the proposed method samples the parameter 
space with the only guide of probability, thus having the 
following advantages: 

• In the fitting procedure, the Bayesian method will 
not get stuck in local minima if its barrier is smaller 
than the error associated with the experimental 
data set. 

• Parameters are obtained as PDFs and, because 
the whole parameter space is sampled, correlations 
between parameters are naturally taken into ac- 
count. Moreover, a natural way to define errors 
based on the PDF of parameters is obtained within 
this method, which following the frequentist ap- 
proximations would be the 68% confidence interval 
around the most probable parameter value, i.e. the 
parameter is inside these limits with a probability 
P=0.68. PDFs may take an arbitrary form, for ex- 
ample indicating that only a superior limit to the 
parameter can be extracted from the experimental 
data. 

• The likelihood (which as we have seen is directly 
related to x 2 ) obtained with this method is also a 
PDF hence revealing the whole complexity of the 
parameter landscape. Model selection is then per- 
formed taking into account all parameter combina- 
tions compatible with the experiment. 

• The presented method is flexible enough to deal 
with low counts experiments where the Poissson 
distribution cannot be approximated by a Gaus- 
sian function, by simply redefining x 2 in the Gibbs 
sampling algorithm. 

This work was supported by the Spanish Ministry of 
Science and Technology (FIS2008-00837) and by the Cat- 
alonia government (2005SGR-00535). 
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